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We show how canonical ensemble expectation values can be extracted from quantum Monte 
Carlo simulations in the grand canonical ensemble. In order to obtain results for all particle sectors, 
a modest number of grand canonical simulations must be performed, each at a different chemical 
potential. From the canonical ensemble results, grand canonical expectation values can be extracted 
as a continuous function of the chemical potential. Results are presented from the application of 
this method to the two-dimensional Hubbard model. 
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I. INTRODUCTION 

The properties of strongly correlated electron systems 
near the Mott insulating phase depend sensitively on the 
dopingji Most simulations of these systems have been 
carried out within the grand canonical ensemble, where 
a convenient formalism exists for the evaluations of finite 
temperature Green's functions and other physical quan- 
tities that can be expressed in terms of themiS In this 
framework, to explore the doping dependence, it is nec- 
essary to carry out simulations at various discrete values 
of the chemical potential {/i Q }, and interpolate between 
them. In this paper we describe a method for optimally 
combining data from these simulations to obtain results 
for a continuous range of (j,. We first evaluate canonical 
ensemble expectation values and partition functions (up 
to an overall constant) from simulations performed in the 
grand canonical ensemble. Using the approach of Ferren- 
berg and Swendsen 3 to combine results from simulations 
performed with different value of the chemical potential, 
we obtain canonical ensemble quantities for a wide range 
of fillings from a modest number of simulations. From 
these we are able to construct grand canonical expecta- 
tion values, enabling us to study a variety of physical 
quantities as a continuous function of the chemical po- 
tential. 

In Section II we present our methodology, and in 
Section III we illustrate it with results for the two- 
dimensional Hubbard model. 



II. METHODOLOGY 

We begin by briefly summarizing the approach to the 
simulation of strongly correlated many-electron systems 
in the grand canonical ensemble set out in Ref. 0- The 
expectation value of a physical observable O is 



(0(M)> 



Tr [0< 



r e -/3(i?- M ;v) 



(1) 



where H is the Hamiltonian, (3 the inverse temperature, 
\x the chemical potential and N the number operator for 
the electrons. In order to perform a numerical simu- 
lation, one must first evaluate the traces over the elec- 
tron degrees of freedom. This is possible if the Hamil- 
tonian is quadratic in the electron creation and annihi- 
lation operators, or can be made so through a Hubbard- 
Stratonovich transformation. To this end we introduce 
a small imaginary-time step, At, by writing — At L. 
The partition function can then be written in the form 

L 



Z{n) = Tr 



= Tr 



,-Ar(H-p,N) 



(2) 



For each time-slice £ — 1 . . . L, we introduce a set of 
Hubbard- Stratonovich variables x(l) such that 



-AtH 



u(x{£)) e- AT Su„ 



f (*(*)) ■ 



(3) 



Here cJ CT and c ia are the creation and annihilation oper- 
ators for electrons at lattice site i with z-component of 
spin <r; /i? • (x(£)) is a single particle Hamiltonian for an 
electron propagating through the external field, x(£); and 
u> (x(£)) is a positive definite weight function. J2x(i) indi- 
cates an integral over continuous Hubbard-Stratonovich 
variables or a sum over discrete ones. Typically, x{£) has 
a component for each spatial lattice site or link. 

The traces in Eq. Q can now be performed, yielding 
an expression of the general form 



(0(m)> 



J2x p{x,/j>)o(x,h) 



(4) 



Here x stands for the totality of Hubbard-Stratonovich 
variables on all time slices, 

L 

p(x, M ) = £ T (x, M ) D l (x, (i) J] uj (x(£)) (5) 

and the determinants for spin up and down electrons are 
given by 

D a (x, (i) = Det [/ + e " " A a (x)] , (6) 
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where / is the unit matrix and 

A <7 (x) = e- ATh °M L »...e- ATh °^ 1 ». 



(7) 



The quantity 0(x, p) in Eq. (0J can generally be ex- 
pressed in terms of finite temperature Green's functions 
for a single electron propagating in the background field 
provided by the Hubbard- Stratonovich variables, x. For 
example, 



Tr 



H cr °j a' c 



D^(x,p) Di(x,fi) 8 at(T > 

X (i+e^A^x)) ■ ( 8 ) 



For models with particle- hole symmetry, such as the Hub- 
bard model at half-filling, the product of the electron 
determinants is positive, and one can use importance 
sampling techniques to generate a sequence of Hubbard- 
Stratonovich configurations with the probability distri- 
bution 



P(x,p) 



p(x,p,) 



(9) 



The average value of 0(x,p) in these configurations is 
then an estimator for (O(p)}. Details of an algorithm for 
efficiently generating configurations are given in Ref. 0. 

For systems which do not have particle-hole symme- 
try, such as the Hubbard model away from half-filling, 
the product of electron determinants will in general not 
be positive definite. In such cases, one can generate 
Hubbard-Stratonovich fields using the probability distri- 
bution 



P u {x,p) 



\p{x,p)\ 



Ex- W^)\ 

It is then necessary to move the sign of p(x,p), 

p(x,p) 



S(x,p) 



±1 



\p(x,p)\ 

into the measurements yielding 

(O(p)) = Pfan)0(x,n) 

X 

J2x p \\(. x ^) Q(x,p) S{x,p) 
E^ p u(x,p) S(x,p) 

The expectation value of the sign can be written 
(S(p)) = E P\\{x,n) S(x,n) 

X 

J2 X p{x,n) = z{p) 
T, x > \p{x',p)\ z n (n)' 



(10) 



(ii) 



(12) 



(13) 



Here Z{p) is the partition function of the physical sys- 
tem, and Z\ | (p) that of a fictitious one in which the sign 
of the product of determinants, S{x,p), is ignored. 



To obtain information about the canonical ensemble 
from grand canonical simulations, we note that so long as 
the electron number operator commutes with the Hamil- 
tonian, and the Hubbard-Stratonovich variables are cho- 
sen so that Eq. © holds, then the product of electron 
determinants has an expansion in the fugacity of the form 



E 

JV 



Z N (x) e 



f3^N 



(14) 



Once we have gone to the computational expense of per- 
form an LDU decomposition of A a (x), Eq. J7|), which 
we must do each time we make a measurement^ it is 
straightforward to evaluate the left-hand side of Eq. (|14|) 
for a number of different values of p. Eq. (|14fl then yields 
a set of linear equations that can be solved for the Zn{x). 
At moderate to low temperatures, only a limited subset 
of the Zn(x) will make a significant contribution to the 
product of determinants, so the system of equations to be 
solved is considerably smaller than the number of spatial 
lattice points. Since the canonical partition function for 
the sector with electron number N is given by 



P\\{x,p) Z N {x) 
^ \D ] {x,p)D i {x,p)\ 



J N 



= Z N (p), (15) 



where 



Zn — ^ z N {x) n w (*(*)) > 



(16) 



we can evaluate Zn using an ensemble of Hubbard- 
Stratonovich fields generated with the probability dis- 
tribution P|| (x, p). 

If the operator O is defined on a single imaginary-time 
slice, or if it does not change the electron number from 
time slice to time slice, then we can also write 

0(x,p) D r (x,p) D L (x,p) = £ O n (x) e^ N , (17) 

and we can obtain a set of linear equations for the On(x) 
by evaluating the left-hand side of Eq. (|17fl for different 
values of p. In this case 



E 



o 



N 



P\\{x,p) On{x) _ 

\Df(x,p) D L (x,p)\ Z n (p) 



O n (p). (18) 



Finally, the expectation value of the operator O in the 
canonical ensemble sector with electron number N is 



(O) 



N 



N 



Ojv 
Zn 



Also, once the Zn and On are in hand, 

E N N e^ N 



(O) 



(19) 



(20) 
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gives the grand canonical expectation values as continu- 
ous functions of \x. 

From simulations at a single value of /i one only expects 
to be able to make accurate determinations of the Zn 
and On for N in the vicinity of (N). We must therefore 
perform a set of simulations with chemical potentials /j. a , 
sufficiently spaced to cover the range of N relevant to the 
problem of interest. As indicated in Eqs. i|15|) and (|18|l . 
the outputs of our simulations are Zn{^ u ) = Zn /Z^^a) 
and On{^o) = On /Z u (/i a ), rather than Zn and On- 
We can combine results from simulations with different 
values of the chemical potential by writing 

Zn = c Q Z N (fJ. a ) Z n (/j, a ) (21) 

a 

O n = J2 d a 5 N (n a ) Z n (n a ), (22) 

a 

c a = J2 d a = 1. (23) 



with 



Following Ferrenberg and Swendsen, we choose the c a 
and d a to minimize the variance of Zn and On subject 
to the constraints of Eq. I|23|l . A short calculation yields 



E 7 VI-Vm-v) a N^i)Y 



(24) 



where <r N (p a ) is the variance of Zjv(/x q ), which we de- 
termine from the simulation. Of course, a corresponding 
result holds for the d a with a N (p a ) replaced by the vari- 
ance of the Oiv(jUa)- 

The constants Z H {fi a ) can be determined up to an over- 
all normalization by iteratively solving the equation 



Z(Ho 



E 

N 



Zn e 



0H*N 



(25) 



with the Z N given by Eqs. (O and The (S^a)) 

are measured directly in the simulations. 

It is also possible to obtain i?n(/x), and therefore 
{S(p)), as a continuous function of [i. We simply note 
that 



(26) 



A' 



where 



ZnM=Y, Zn(x) S(x,fi)Y[ u(x(£)). (27) 



Note that once the Zn(x) are known for any configura- 
tion, we can determine S(x, fi) for any value of /i from the 
right-hand side of Eq. (|14l) . A simulation performed at 
the chemical potential, will, of course, only determine 
the ratio Zn Z H {fi a ). However, we can combine re- 
sults from simulations performed at different values of 
the chemical potential just as for the Zn- 
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FIG. 1: Free energy difference between half-filled system and 
system with chemical potential fi. Statistical errors are neg- 
ligible on this scale. 



III. NUMERICAL RESULTS 

We illustrate the methodology outlined in the last sec- 
tion with results for the two-dimensional Hubbard model. 
The Hamiltonian is 



H 



<ij>,a 



~ C ja C i<7' 



-^E 



Om-DOu-s)- 

(28) 

Here c icr and c ia are the creation and annihilation oper- 
ators for electrons with z-componcnt of spin a at lattice 
site i, and ni a — c i ia c ia . The sum < ij > is over all 
pairs of nearest neighbor lattice sites, t is the hopping 
parameter, and U the Coulomb coupling constant. 

The Coulomb term is reduced to quadratic form in 
the electron creation and annihilation operators using 
Hirsch's discrete Hubbard- Stratonovich transformation 4 



e -Ari7(n n -i)(n u -i) = 



(29) 



For /i = 0, which corresponds to half- filling, particle- 
hole symmetry implies that D^(x, 0) D^(x, 0) is always 
positive, 4 so S(x, 0) = 1, and there is no sign prob- 
lem. It is therefore convenient to adopt the normalization 
Z n (0) = Z(0) = 1 in solving Eq. for Z n (fj, a ). Thus, 
we are in fact able to use Eqs. (|25|l and H26(l to determine 
Z(fj,)/Z(0) and Z n (fj,)/Z u (0) respectively. 

All of the results we present here are on a 4 x 4 lattice 
with t = 1 and U = 4. The number of time slices, L, 
is chosen so that At = 1/8. Except where otherwise 
noted, simulations were performed at ji = —1.025 and 
fi = —0.6 for each temperature. At (3 = 8 we performed 
additional runs at both fi = —0.9625 and [i — —0.9, while 
at other temperatures either /i = —0.9625 or /i = —0.90 
was used. For runs at fi = —0.6, 100,000 Monte Carlo 
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FIG. 2: Density of the system for several different values 
of p. As errors depend on (i, errorbars are shown at several 
points along the curves here and in subsequent figures. Also 
shown is the zero-temperature result, calculated using exact 
diagonalization. 
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FIG. 3: Charge compressibility of the system, obtained by 
the analytical differentiation of n. 



sweeps with 10,000 warmup sweeps were performed. For 
all other runs, 400,000 Monte Carlo sweeps with 10,000 
warmup sweeps were performed. For all simulations, non- 
local moves, as suggested in Ref. y, were used to assure 
ergodicity. To invert Eqs. |Q and (|T7jl . the right-hand 
sides are measured at the set of chemical potentials, 



(30) 



where [i a is the chemical potential used in the simu- 
lation, i = — 7...0...7 and <5/x = 0.02. After inver- 
sion and averaging over configurations, particle sectors 
where Z^e^ >laN /Z(/j, a ) < 10~ 4 are dropped to prevent 
the spread of roundoff error from the inversion. The jack- 
knife method was used for error analysis. It should be 
noted that, after analysis, results at different values of /j, 
are not statistically independent. 

In Fig. n we plot the free energy difference, 



FQjl) - F(0) 



(31) 



as a function of /i, for two values of (3. In Fig. [21 we plot 
the density defined by n(fj,) = (N)/V. Here V is the 
number of spatial lattice points and (N) is calculated 
using the standard thermodynamic identity, 



(N) 



dF 



Y,n n Z n e 
J2n z n e 



0»N 



(32) 



As the temperature is lowered, the transition between the 
half-filled state (n — 1.0) and the 6-hole state (n = 0.625) 
becomes sharper. In particular, at zero-temperature the 
density decreases in a series of jumps, due to the discrete- 
ness of the finite-size spectrum. 

Within our framework it is also straightforward to 
compute the compressibility of the system, k = dn/dfj,, 
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FIG. 4: (a) The antiferromagnetic structure factor, S zz (ji, it) . 
(b) The d-wave pair field correlation function. 
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FIG. 5: The expectation value of the sign. Calculated from 
runs with 100,000 Monte Carlo sweeps at /i = —1.5, \i = 
-1.025, and fi = -0.6. 

by differentiation of Eq. i(32() . Note that the differen- 
tiation can be performed analytically. Figure [3] shows 
the compressibility as a function of /i for different values 
of f3. As the temperature is lowered, the compressibil- 
ity develops a peak around fi = —1.0 that is likely to be 
the signature of the low-temperature divergence expected 
from the metal-insulator transition^ 

As previously mentioned, within our numerical 
scheme it is possible to calculate observables that are 
not diagonal in the particle number. Figure Ufa) shows 
the antiferromagnetic structure factor. This is given by 

S^7r,7r)=i£(-l) i+ ''S?S;, (33) 
y 

where Sf = \c\ a a z afi ap is the standard spin operator. 
The plot of this quantity vs. /i clearly indicates that 
the antiferromagnetic correlations present at half-filling 
are sharply suppressed upon doping. A similar plot of 
the equal-time d-wave pair field correlation function is 
shown in Fig.^Jb). Here the rf-wave pair field correlation 
function is given by 



where Aj = \ J2s(~ 1) <Sc 1t C '!+<U creates two electrons in 
a d-state. Here, S sums over the four near-neighbor sites 
of / and (— l) s gives the sign alternation characteristic 
of a d-wave pairing amplitude. The enhancement of p c i 
toward /i — 0.0 is a finite-size effect due to a strong 
antiferromagnetic response in the nearest-neighbor terms 
in Eq. (jMjl . 

Finally, we show the expectation value of the sign in 
Fig. [SJ Here the sign is calculated as a continuous func- 
tion of /i using Eqs. 1)26(1 and ((27(1 . Note that the sign is 
small in the fi = — 1.0 region where the density is chang- 
ing rapidly and electron correlations are believed to be 
important. The statistical fluctuations of the other ob- 
servables grow as the sign decreases. 

IV. CONCLUSION 

In this paper we have presented a new method for ex- 
tracting canonical ensemble results from grand canonical 
ensemble quantum Monte Carlo simulations. As canoni- 
cal information is only extracted from sectors whose par- 
ticle number is close to the average number of particles 
in the simulation, simulations must be performed at sev- 
eral different chemical potentials to obtain results for a 
range of particle number sectors. These separate simula- 
tions can then be combined to obtain a complete picture 
of the different canonical ensembles with lower statistical 
fluctuations than any of the simulations taken individu- 
ally. Once the canonical results are obtained, they can be 
combined to give grand canonical results as a continuous 
function of the chemical potential. 

In this work we have presented results for the two- 
dimensional Hubbard model on a 4 x 4 lattice with 
Coloumb interactions of moderate strength, but the 
method is applicable to any quantum mechanical prob- 
lem, simulated in the grand canonical ensemble, for which 
particle number is conserved. 
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